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A method for calculating the Kohn-Sham exchange-correlation potential, Vxc( r ), from a given 
electronic wavefunction is devised and implemented. It requires on input one- and two-electron 
density matrices and involves construction of the generalized Fock matrix. The method is free from 
numerical limitations and basis-set artifacts of conventional schemes for constructing vxc(r) in which 
the potential is recovered from a given electron density, and is simpler than various many-body 
techniques. The chief significance of this development is that it allows one to directly probe the 
functional derivative of the true exchange-correlation energy functional and to rigorously test and 
improve various density-functional approximations. 
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The Kohn-Sham density-functional theory [1] is the 
most widely used method for electronic structure cal¬ 
culations of molecules and solids. In this method, the 
ground-state energy of a system is treated as a func¬ 
tional of the electron density p(r) and then partitioned in 
such a way that only one term, the exchange-correlation 
energy Exc[p]i remains unknown. Application of the 
variational principle to the total energy functional leads 
to a one-electron Schrodinger equation with an effective 
Hamiltonian that includes the system’s electrostatic poten¬ 
tial and the exchange-correlation potential, i’xc([p]; r ) = 
5Exc[p]/3p( r )- While the exact Exc[p\ can be writ¬ 
ten only in implicit form [2], its functional derivative 
■uxc([p];r) can in principle be computed and visualized 
as a function of r for any particular non-interacting v- 
representable density. High-quality Kohn-Sham poten¬ 
tials are used for testing density-functional approxima¬ 
tions [3], accurate description of electronic excitations [4], 
and other purposes. 

Most existing methods for generating exact exchange- 
correlation potentials fit the function t>xc( r ) to a given 
ground-state p( r) via the Kohn-Sham equations either 
by iterative updates [5-8] or through some constrained 
optimization [9-11]. The target densities are usually ob¬ 
tained from ab initio wavefunctions which are themselves 
discarded. Because small changes in p{ r) can induce large 
changes in r>xc( r ) [12], potential-reconstruction meth¬ 
ods that use only p( r) as input suffer from numerical 
instabilities. Moreover, electron densities generated using 
ubiquitous Gaussian basis sets correspond to exchange- 
correlation potentials that wildly oscillate and diverge [13- 
16], a result that is formally correct but unwanted. Kohn- 
Sham potentials can be also constructed by many-body 
methods [17-21], but these techniques are quite elaborate 
and often require solving an integral equation for , uxc( r )i 
which is a challenge by itself. 


Here, we propose a radically different method for 
computing the exchange-correlation potential of a given 
many-electron system, which avoids the above difficul¬ 
ties. In this method, the functional derivative of the 
exact Exc[p) is obtained directly from the system’s elec¬ 
tronic wavefunction. The approach represents a non¬ 
trivial generalization of our technique for constructing 
Kohn-Sham potentials corresponding to Hartree-Fock 
(HF) electron densities [22, 23] and is conceptually re¬ 
lated to the wavefunction-based analysis of Kohn-Sham 
potentials developed by Baerends and co-workers [24-28]. 

The basic idea of our approach is to derive two ex¬ 
pressions for the local electron energy balance, one of 
which originates from the Kohn-Sham equations, the 
other from the Schrodinger equation. When one expres¬ 
sion is subtracted from the other under the assumption 
that the Kohn-Sham and wavefunction-based densities 
are equal, the system’s electrostatic potentials cancel out 
and the difference gives an explicit formula for uxc( r )- 
For simplicity, the treatment presented in this Letter is 
restricted to electronic singlet ground states described 
with closed-shell Kohn-Sham determinants, and assumes 
that all basis functions and orbitals are real (although 
the notation for complex conjugate is retained). 

Accomplishing the first part of this plan is easy. In the 
Kohn-Sham scheme, the ground-state density of a singlet 
IV-electron system is obtained as p KS (r) = JT m\(j)i(r)\ 2 , 
where rq = 0 or 2 are occupation numbers of the corre¬ 
sponding Kohn -Sham orbitals (N = W ; n,). The orbitals 
are obtained by solving the equation 



+ u(r) + 'CH S (r) + «xc( r ) 




£i0i(r), (1) 


where u(r) is the electrostatic potential of the nuclei 
and Uy S (r) is the electrostatic potential of p KS (r). If 
we multiply Eq. (1) by r), sum over i , and divide 
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through by p KS (r), we obtain 


and 


-r KS fri 

^KS( r ) + w ( r ) + ^H S ( r ) + ^xc(r) = e KS (r), ( 2 ) 

where r]j s ( r) = —(1/2) J2i n i^t (r)V 2 (/i(r) is the Kohn- 
Sham kinetic energy density and 

" KS(r) = piy E n M<t>i{*)\ 2 (3) 

is the average local Kohn-Sham orbital energy [29] . 

The second part of the plan is to reduce the IV-electron 
Schrodinger equation to a local energy balance expression 
analogous to Eq. (2). There is more than one way to do 
this. Holas and March [30] had considered a contracted 
Schrodinger equation for this purpose, but their proposal 
led to a complicated integral equation for uxc( r ) involv¬ 
ing the three-particle reduced density matrix (RDM). 
The Baerends group [24-28] used an expression involving 
(N — l)-electron conditional amplitudes. The method 
we propose here is motivated by Lowdin’s approach [31] 
to the problem of finding the optimal finite one-electron 
basis set for a configuration interaction (Cl) expansion. 

Suppose we have an iV-electron ground-state wavefunc- 
tion T expressed in terms of orthonormal orbitals (i/'i}. 
Then the total electronic energy may be written as 

E = E lij^jVA^i) \ r \2 IV’iV’fe)) (4) 

ij ikjl 


where h( r) = — (1/2)V 2 + v(r) is the one-electron core 
Hamiltonian, 7 ^ = X/ 0 .(’4 , |dj 0 .aio-|'I') (<r = a,/3 is the spin 
index) are matrix elements of the spin-free one-particle 
RDM, and T ikj i = (1/2) (’I'|a] (T aj CT ,a fc(T /a^ CT |1 , ) are 
matrix elements of the spin-free two-particle RDM. 

Our objective is to turn Eq. (4) into a local energy 
balance equation. We start by minimizing E with respect 
to the functions {ipi}, subject to the constraint ( ipj\ipi ) = 
Sji, while keeping 7y and Ti k ji fixed. The corresponding 
Euler-Lagrange equation is 

(5) 

where A ij are yet undetermined Lagrange multipliers. We 
evaluate the functional derivative in Eq. (5), multiply the 
result by i/i* (r'), sum over j, and obtain 

ft(r)7( r ,r , )+2 j dr ^ = E A ^i( r W( r 0 - 

1 1 ij 

( 6 ) 

where 


r(r,r 2 ;r , ,r , 2 ) = r »kj^i( r )V’*( r 2 )V’j {r')ipl{r' 2 ) ( 8 ) 

ikjl 

are the coordinate representations of the spin-free one- 
and two-particle RDMs, respectively. 

We denote the left-hand side of Eq. ( 6 ) by G(r, r') and 
treat it as the kernel of an integral operator defined by 

G'tpjir) = J G(r, r')ipj(r') dr'. (9) 


Then X,j can be determined from Eqs. ( 6 ) and (9) as 


a a = (10) 


The operator G, known as the generalized Fock operator or 
orbital Lagrangian, arises in various problems of quantum 
chemistry [31-35]. 

For our purposes, we need only the r = r' part of Eq. ( 6 ) 
which after division by p WF (r) = y(r, r) becomes 


(r) 


WF 

—— 1 —- -t- c(r) -|- 

^ WF (r) v ' „wf 


( r ) 


P(r, r 2 ) 

|r-r 2 | 


dr 2 = e w *(r), ( 11 ) 


where r/ VF (r) = —(1/2) [V 2 y(r, r')] r , =r is the interacting 
kinetic energy density, P(r, r 2 ) = T(r, r 2 ; r, r 2 ) is the pair 
function, and 


-WF 


( r ) = 7WF)viE A b-^( r )V'j*(i')- (12) 


One can always write the pair function as 

p (r, r 2 ) = ^p WF ( r ) [p WF (r 2 ) + Pxc (r, r 2 )] , (13) 


which defines Pxc'( r , r 2 )> the exchange-correlation hole 
density. Substituting Eq. (13) into Eq. (11) we obtain 


( r ) 


r WF 

pWF(r) + «(r) + < F (r) + ^ F (r) = e WF (r), (14) 


where u^ F (r) is the electrostatic potential of p WF (r) and 



(15) 


is the Slater exchange-correlation-charge potential [36]. 
Equation (14) is the wavefunction counterpart of Eq. (2). 

Observe that the sum in Eq. (12) does not change if 
we replace every X,j with X* { . This means that e WF (r) 
is determined by the Hermitian (symmetric) part of G. 
If desired, one can define the self-adjoint operator F = 
(G + CP)/2 and solve the Hermitian eigenvalue problem 
Ffi{ r) = A,/,;(r). This optional step allows one to cast 
Eq. (12) as 


7 ( 1 -, r') = ^ 7 ij V’i(r)^*(r') 
ij 


1 


n WF ( r 


E A *i^( r )i 2 > 


( 7 ) 


-WF / \ 

e (r) 


(16) 
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which is formally analogous to Eq. (3). The quantity 
e WF (r) as given by Eq. (16) was introduced by us earlier 
under the name of “average local electron energy” [37]. 

Now let us subtract Eq. (14) from Eq. (2), sub¬ 
stitute the identity tl = t — 'S7 2 p/‘i for t^ s 
and for t/ vf with r KS = ( 1 / 2 ) n i|V 0 i | 2 and 
r WF (r) = (1/2) [V r 'V r 7 (r, r / )] r , =r , and apply the con¬ 
dition p KS (r) = p WF (r). This yields the central equation 
of this work: 


vxc(r) = TT’M + < 


3 (r) - e WF (r) 


r WF 


(r) 


rKS, 


r) 


pWF ( r ) p KS ( r )' 

(17) 

Since r KS and e KS are initially unknown, Eq. (17) must 
be solved iteratively in conjunction with the Kohn -Sham 
equations. The transition from tl to r is not strictly 
necessary but beneficial for numerical calculations because 
r does not diverge at the nuclei as does 77,. 

Note that as r —► 00 , the term v^ F vanishes, but the 

? ks ^ks/„ks and 


other ingredients remain nonzero: e 


7p 


— t ks //9 ks approach chomo [38], while e WF , t/ VF /p WF , 
and — r WF /p WF approach — I min [37], where I m i n is the 
first ionization energy of the system as determined by the 
extended Koopmans theorem [39]. To ensure that vxc( r ) 
as given by Eq. (17) properly vanishes at infinity, we shift 
all current values of e; in each Kohn-Sham iteration to 
satisfy the condition 


eHOMO = — Tmin, (18) 

which also imparts p KS (r) with proper asymptotic decay. 

The proposed algorithm is as follows. 

1. Obtain a wavefunction for the system of interest. 
Calculate p WF , r WF , v^ F , e WF , and / m i n - 

2. Generate an initial guess for the occupied Kohn- 
Sham orbitals {4>i} and their eigenvalues {ei}. 

3. Using the current guess for {</>,} and shifted {e,}> 
construct the potential vxc by Eq. (17). 

4. Solve the Kohn-Sham equations using the current 
vxc and the same basis as in step 1. This gives new 
sets {(/>,;} and {ej. 

5. Return to step 3 and iterate until the potential vxc 
is self-consistent. 

The method was implemented in the GAUSSIAN 09 
suite of programs [40] , which already contains subroutines 
for constructing the generalized Fock matrix as part of 
the multiconfigurational self-consistent held (MCSCF) 
module. The values of J m j n were computed as in Ref. 34, 
while p WF and r WF were assembled from natural orbitals. 
Any reasonable density-functional approximation may 
be used to generate an initial guess for {</>*} and {&;}■ 
The potential was considered converged when all Kohn- 
Sham density matrix elements from consecutive iterations 
differed by less than 10 -10 in the root-mean-square sense. 



FIG. 1. Exchange-correlation and correlation (inset) potentials 
for the He atom calculated from FCI wavefunctions using 
various basis sets. 



FIG. 2. Exchange-correlation potentials for the Ne and Be 
atoms calculated from compact CASSCF wavefunctions using 
various basis sets. 

The method works best with basis sets that are not heavily 
contracted in the core region. 

An added benefit of generating uxc( r ) from a wave- 
function is that one can readily obtain the corresponding 
exchange-correlation energy, which is inaccessible 

when only the electron density is known. We computed 
this energy as E^, = E^ F + T c , where E^ F is the 
ab initio exchange-correlation energy defined as E^ F = 
(1/2) J p WF (r)ug VF (r) dr and T c = T — T s is the differ¬ 
ence between the ab initio and Kohn-Sham total kinetic 
energies, evaluated analytically. Also of interest is the in¬ 
tegrated density difference, A p = J |p KS (r) — p WF (r)| dr, 
evaluated for the self-consistent vxc( r )- Because the con¬ 
dition p KS (r) = p WF (r) is imposed in our approach only 
in the derivation of Eq. (17), A p strictly vanishes only 
in the basis-set limit. Insistence on reproducing p WF (r) 
exactly in Gaussian basis sets would be misplaced because 
(i) it brings out unwanted oscillations and divergences of 
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TABLE I. Characteristics of selected wavefunctions and the corresponding Kohn-Sham effective potentials (in atomic units). 


System 

Wavefunction 

Etot 

-^min 

T s 

ii 

4 

1 

4 

to J 

rpKS 

^XC 

A p 

He 

FCI/cc-pVTZ 

-2.900 232 

0.9013 

2.8571 

0.0435 

-1.0550 

0.002 51 


FCI/cc-pVQZ 

-2.902411 

0.9026 

2.8652 

0.0370 

-1.0645 

0.000 65 


FCI/cc-pV5Z 

-2.903152 

0.9032 

2.8668 

0.0364 

-1.0662 

0.00013 


Exact a 

-2.903 724 

0.9037 

2.8671 

0.0366 

-1.0667 


Be 

CAS(2,4)/cc-pCVDZ 

-14.61545 

0.3485 

14.4901 

0.1333 

-2.6146 

0.01729 


CAS(2,4)/cc-pCVTZ 

-14.616 53 

0.3489 

14.5538 

0.0619 

-2.6866 

0.004 93 


CAS(2,4)/cc-pCVQZ 

-14.616 77 

0.3490 

14.5910 

0.0258 

-2.7232 

0.005 47 


FCI/u-cc-pCVTZ 

-14.663 70 

0.3421 

14.5956 

0.0654 

-2.7715 

0.00215 


Exact a 

-14.66736 

0.3426 

14.5942 

0.0732 

-2.7701 



a Accurate estimates from Ref. 41 (He) and Ref. 42 (Be). 


uxc(r) and (ii) the potential that yields a given density 
in a finite basis is not unique anyway [43, 44]. 

To test the method, we computed exchange-correlation 
potentials for the three atoms (He, Be, and Ne) for which 
exact potentials are available in the literature [41, 42] 
using full configuration interaction (FCI) and complete ac¬ 
tive space (CAS) SCF wavefunctions and standard Gaus¬ 
sian basis sets [45]. For He, already the potential extracted 
from the FCI wavefunction in the cc-pVTZ basis set is very 
close to the exact uxc( r ), and the cc-pVQZ and cc-pV5Z 
FCI exchange-correlation potentials are visually indistin¬ 
guishable from the benchmark (Fig. 1 and Table I). Even 
the correlation potential for He, uc(r) = i>xc( r ) —i’H(r)/2, 
which is almost two orders of magnitude smaller than 
vxc(r), is very accurate at the FCI/cc-pV5Z level (Fig. 1). 
For Be, the sequence of potentials from CAS(2,4) wave- 
functions quickly approaches the exact uxc( r ) with in¬ 
creasing basis set size (Fig. 2), as do the corresponding 
T a values (Table I). By contrast, T c and E^c converge 
slowly because they depend not only on uxc( r ) but also 
on the accuracy of the wavefunction through the value of 
T. Potentials for the Ne atom constructed from CAS(8,8) 
wavefunctions also improve rapidly with the size of the 
basis set (Fig. 2). Thus, even compact correlated wave- 
functions can produce accurate Kohn-Sham potentials, 
provided that the basis set is of good quality. 

The method works equally well for molecules. It is 
known that, in molecules, the onset of strong correlation 
induced by bond stretching manifests itself in character¬ 
istic mid-bond peaks of uxc( r ) [27, 46-48]. Using our 
method, we readily reproduced these peaks in a number of 
stretched diatomics exemplified by N 2 (Fig. 3). Exchange- 
correlation potentials for polyatomic molecules can also 
be generated by our method (Fig. 4). 

It is remarkable that Kohn-Sham potentials computed 
from wavefunctions are always well-defined and free from 
spurious features. Conventional methods for extracting 
vxc(r) from densities, when implemented in matrix form, 
would not deliver such unambiguous results because there 
is no one-to-one correspondence between densities and 
potentials in finite basis sets [43]. Furthermore, when 



2 foj 

FIG. 3. Exchange-correlation potentials for the N 2 molecule 
obtained from HF and valence CASSCF wavefunctions at the 
experimental equilibrium bond length and at 2 R e . 



2 hj 


FIG. 4. Exchange-correlation potentials for HCN obtained 
from HF and valence CASSCF wavefunctions at the experi¬ 
mental equilibrium geometry and with A(HC) = 2_R e (HC). 


density-to-potential mapping techniques are rigorously 
applied to electron densities generated in Gaussian basis 
sets, one obtains unphysical potentials [13-16]. Neither 
of these complications affects our approach. 
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In conclusion, we have developed a practical method for 
folding a many-electron wavefunction into the correspond¬ 
ing exchange-correlation potential. The key ingredient 
of our approach is the generalized Fock matrix which is 
commonly available in ab initio codes as a by-product 
of computing MCSCF wavefunctions, nuclear gradients, 
and first-order properties. The method possesses sev¬ 
eral advantages over existing techniques for constructing 
exchange-correlation potentials: it delivers cxc( r ) i n a 
simple analytic form, avoids the ambiguity of associating 
a given electron density with a Kohn-Sham potential in 
a finite basis set, is stable with respect to changes in 
basis sets, convergence thresholds and other details of the 
calculation, and produces potentials without oscillations 
and divergences when using Gaussian basis sets. Further 
exploration of the capabilities of our approach is under 
way. 
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